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Numerical  Determination  of  Minimum  Mass  Structures 
with  Specified  Natural  Frequencies^ 

by 

A.  MIELE^,  A.  MANGIAVACCHI^,  B.P.  MOHANTy'*  , and  A.K.  WU^ 

Abstract.  The  problem  of  the  axial  vibration  of  a cantilever 
beam  is  investigated  both  analytically  and  numerically.  The 
mass  distribution  that  minimizes  the  total  mass  for  a given 
value  of  the  frequency  parameter  3 is  determined  using  both 
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the  sequential  ordinary  gradient-restoration  algorithm  (SOGRA) 
and  the  modified  quasilinearization  algorithm  (MQA) . Concern- 
ing the  minimum  value  of  the  mass,  SOGRA  leads  to  a solution 
precise  to  at  least  4 significant  digits  and  MQA  leads  to  a 
solution  precise  to  at  least  6 significant  digits. 

Comparison  of  the  optimal  beam  (a  variable-section  beam) 
with  a reference  beam  (a  constant-section  beam)  shows  that  the 
weight  reduction  depends  strongly  on  the  frequency  parameter  8. 
This  weight  reduction  is  negligible  for  3-*-0,  is  11.3%  for 
3=1,  is  55.3%  for  3 = 1.4,  and  approaches  100%  for  ^-*-'rr/2. 

Key  Words.  Structural  optimization,  dynamic  optimization, 
axial  vibrations,  frequency  constraint,  fundamental  frequency 
constraint,  optimal  structures,  cantilever  beams,  bars,  sequen- 
tial gradient-restoration  algorithm,  modified  quasilinearization 
algorithm,  numerical  methods,  computing  methods. 
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1.  Introduction 

Under  the  sponsorship  of  the  Office  of  Scientific  Re- 
search, Office  of  Aerospace  Research,  United  States  Air  Force, 
the  Aero-Astronautics  Group  of  Rice  University  has  developed 
in  recent  years  several  algorithms  dealing  with  the  numerical 
solution  of  optimal  control  problems  on  a digital  computer. 

Noteworthy  among  those  algorithms  are  the  sequential 
ordinary  gradient-restoration  algorithm  (SOGRA,  Refs.  3-4)  and 
the  modified  quasilinearization  algorithm  (MQA,  Refs.  5-6) . 

For  a survey  of  these  algorithms,  see  Ref.  7. 

The  algorithms  discussed  in  Refs.  3-7  are  general,  all- 
purpose algorithms.  As  far  as  aerospace  engineering  is  con- 
cerned, these  algorithms  can  be  applied  to  a variety  of  fields, 
for  example,  optimum  trajectories,  optimum  aerodynamic  shapes, 
and  optimum  structures.  With  reference  to  optimum  structures, 
two  categories  of  problems  can  be  identified:  (i)  static  op- 
timization problems,  for  example,  minimization  of  the  weight 
of  a structure  for  a given  load  distribution;  and  (ii)  dynamic 
optimization  problems,  for  example,  minimization  of  the  weight 
of  a structure  for  given  frequency  constraints. 

This  paper  is  the  first  of  a series  dealing  with  dynamic 
optimization  problems  and  is  being  written  at  the  suggestion 
of  Dr.  V.B.  Venkayya,  Flight  Dynamics  Laboratory,  Wright 
Patterson  AFB,  Ohio.  The  questions  posed  by  Dr.  Venkayya  were 
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as  follows:  can  SOGRA  and  MQA  be  usefully  applied  to  the 
solution  of  dynamic  optimization  problems?  How  reliable  and 
precise  are  the  solutions  obtained  with  these  algorithms  under 
realistic  engineering  conditions? 

At  Dr.  Venkayya's  urging,  the  problem  of  the  axial 
vibration  of  a cantilever  beam  was  investigated,  with  the 
following  objective  in  mind;  to  find  the  mass  distribution 
that  minimizes  the  total  mass  for  a prescribed  fundamental 
frequency.  This  problem  was  studied  first,  because  a closed- 
form  solution,  due  to  Turner  (Ref.  8),  is  available.  Hence, 
it  constitutes  an  ideally  suited  test  problem  for  checking 
the  accuracy  and  reliability  of  SOGRA  and  MQA  in  the  dynamic 
optimization  of  structures. 

Outline.  This  paper  is  divided  in  three  parts:  physi- 
cal problem  (Sections  2-4),  optimization  problem  (Sections 
5-6) , and  numerical  results  (Sections  7-9) . The  first  part 
contains  a description  of  the  physical  problem,  both  in  di- 
mensional form  (Section  2)  and  in  dimensionless  form  (Section 
3);  a reference  solution,  that  pertaining  to  a constant-section 
beam,  is  given  in  Section  4.  The  second  part  contains  a des- 
cription of  the  mathematical  problem  (Section  5) ; Section  6 
presents  several  alternative  formulations  and  an  essential 
simplification;  the  elimination  of  a redundant  constraint. 

The  third  part  contains  the  description  of  the  algorithms 
(Section  7),  the  experimental  conditions  (Section  8),  the 
numerical  results  (Section  9) , and  the  conclusions  (Section  10). 
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Notation.  Throughout  the  paper , the  following  notation 
is  employed. 

Physical  System 

2 

A Cross-sectional  area,  ft 

E Modulus  of  elasticity,  lb  ft  ^ 

F Elastic  force,  lb 

e 

Inertia  force,  lb 
L Length  of  the  beam,  ft 

-2  2 

M Mass  per  unit  length,  lb  ft  sec 

-1  2 

M Tip  mass,  lb  ft  sec 

o 

-1  2 

M*  Total  mass  of  the  beam,  lb  ft  sec 

t Time,  sec 

U Axial  displacement,  ft 

X Axial  coordinate,  ft 

y(x)  Displacement  function,  ft 

Z(t)  Temporal  function  , dimensionless 

e Normal  strain  , dimensionless 

-4  2 

p Density,  lb  ft  sec 

-2 

c Normal  stress,  lb  ft 

0)  Natural  frequency,  sec  ^ 

Subscripts  for  the  Physical  System 

X Derivative  with  respect  to  the  axial  coordinate 

t Derivative  with  respect  to  time 
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Normalized  System 

I Total  mass  of  the  beam,  I = M*/M 

* o I 

m Mass  per  unit  length,  m = ML/M^ 

t Axial  coordinate,  t = x/L 

u First  derivative  of  the  displacement,  u = y I 

w Second  derivative  of  the  displacement,  w=y 

y Axial  displacement,  y = Y(x)/y(L) 

z Auxiliary  variable,  z = mu  i 

6 Frequency  parameter,  B = wL/{p/E) 

Superscript  for  the  Normalized  System 

Derivative  with  respect  to  the  axial  coordinate  t 
(for  example,  y = dy/dt) 
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2.  Physical  System 

We  consider  the  problem  of  the  axial  vibration  of  a can- 
tilever beam  or  bar,  that  is,  a beam  having  a fixed  end  and  a 
free  end.  Attached  to  the  free  end  is  the  concentrated  mass 
. The  beam  has  variable  cross  section;  hence,  it  has  va- 
riable mass  per  unit  length  M.  The  total  mass  of  the  beam  is 
given  by 


Mdx  , 


where  x denotes  an  axial  coordinate,  measured  from  the  fixed 
end  of  the  beam,  and  L is  the  length  of  the  beam.  Note  that 
X = 0 at  the  fixed  end  and  x = L at  the  free  end. 

Assume  that  the  system  composed  by  the  masses  M and  M* 
is  excited  by  some  external  axial  disturbance.  Then,  the  beam 
vibrates  axially  under  the  combined  effect  of  the  elastic  for- 
ces and  the  inertia  forces.  These  forces  are  in  equilibrium. 

Equilibrium  Equation.  Let  denote  the  elastic  force  at 

station  x;  let  F denote  the  elastic  force  at  station  x=x  + dx; 
e 

let  F^  - F^  denote  the  net  elastic  force  acting  on  the  beam 
element  having  length  dx;  and  let  dF^  denote  the  inertia  force 
associated  with  the  mass  element  Mdx.  The  equilibrium  condi- 
tion between  elastic  forces  and  inertia  forces  requires  that 


F - F = dF . , 

e e 1 


Note  that  M = M (x) . 
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where 
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Fg  = aA,  (3-1) 

Fg = oA + [ 3 (aA)/8xldx  , (3-2) 

dF^  = MU^^dx  . (3-3) 

In  the  above  equations,  a denotes  the  normal  stress,  U the  axial 
displacement,  the  axial  velocity,  the  axial  acceleration, 

and  t is  the  time.  From  (2)- (3),  the  equilibrium  condition  takes 
the  form 

3 (aA)/3x  = . (4) 

The  left-hand  side  of  Eq.  (4)  can  be  expressed  in  terms 
of  M and  U if  two  observations  are  made.  First,  the  mass  per 
unit  length  M,  the  cross-sectional  area  A,  and  the  density  of 
the  material  p satisfy  the  relation 

M=  pA  . (5) 

Second,  on  account  of  Hooke's  law,  the  normal  stress  a,  the 

normal  strain  e = U^,  and  the  modulus  of  elasticity  E satisfy 
0 

the  relation 


0 = Ee  = EUj^  . (6) 

Therefore,  upon  combining  Eqs.  (4)-(6),  the  equilibrium  condition 

^Note  that  A = A(x). 

0 

In  Eq.  (6),  a one-dimensional  state  of  stress  is  assumed. 
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takes  the  form 

(E/p)  0/3x)  ( MU^)  - MU^^=  0 . (7) 

Boundary  Conditions.  The  partial  differential  equation 
(7)  must  be  solved  subject  to  appropriate  boundary  conditions. 

(i)  At  the  fixed  end  of  the  beam,  the  displacement  must 
vanish.  Therefore, 


U = 0 at  X = 0 . (8) 

(ii)  At  the  free  end  of  the  beam,  the  elastic  force  act- 
ing on  the  cross  section  A must  balance  the  inertia  force 
associated  with  the  concentrated  mass  . Clearly, 

Fg  + = 0 at  X = L , (9) 

where 

Fg  = aA=  (E/p)MU^  , (10) 


F . 
1 


= M U.  . 
o tt 


Therefore,  upon  combining  Eqs.  (9)-(ll),  we  see  that 

(E/p)MU  + M = 0 at  x = L . 

X o tt 


(11) 


(12) 


(iii)  The  spatial  boundary  conditions  (8)  and  (12)  must 
be  completed  by  appropriate  temporal  boundary  conditions. 


10 


AAR-138 


However,  these  temporal  boundary  conditions  are  not  specified 
here:  we  are  only  interested  in  the  spatial  behavior  of  the 
system. 

Separation  of  Variables.  We  assume  that  the  nature  of 
the  temporal  boundary  conditions  is  such  that  the  technique  of 
separation  of  variables  is  applicable  to  the  system  composed 
of  Eqs.  (7),  (8),  (12).  Therefore,  we  write  the  unknown  function 
U(x,t)  in  the  form 

U(x,t)  = Y(x)  Z (t)  . (13) 

Upon  substitution,  we  see  that  the  previous  differential  system 
splits  into  a spatial  system  and  a temporal  system,  each  des- 
cribed by  ordinary  differential  equations. 

(i)  For  the  spatial  system,  Eqs.  (7),  (8),  (12)  can  be  re- 
written as 

(E/p)  (d/dx)  (MY^)  + = 0 , (14) 

Y = 0 at  X = 0 , (15) 

(E/p)MY  - (o^M  Y=0  at  X = L , (16) 

X o 

where  u>  denotes  any  of  the  natural  frequencies  of  the  system 

9 

under  consideration.  Clearly,  the  solution  Y(x)  of  (14) -(16) 
depends  on  the  mass  distribution  M(x)  as  well  as  on  several 
constants  characteristic  of  the  system  (L,  E,  p ,aj  ,M^)  . 

9 2 

Note  that  w is  the  separation  constant,  that  is,  the  constant 
which  allows  the  separation  of  the  spatial  system  from  the 
temporal  system. 
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(ii)  For  the  temporal  system,  Eq.  (7)  implies  that 

= 0 . (17) 

Of  course,  this  equation  must  be  completed  by  appropriate  and 
consistent  temporal  boundary  conditions. 

(iii)  In  the  remainder  of  this  paper,  we  focus  our  atten- 
tion on  the  solution  of  the  spatial  system  (14)-(16),  For 
problems  where  the  function  M(x)  is  prescribed  a priori  (for 
instance,  the  study  of  the  axial  vibration  of  a constant  section 
beam),  the  system  (14) -(16)  is  linear.  For  problems  where  the 
function  M(x)  is  to  be  determined  (for  instance,  the  optimiza- 
tion problem  described  later  on  in  this  paper) , the  system 

(14) -(16)  is  nonlinear.  In  either  case,  the  system  (14) -(16) 
is  homogeneous  in  the  unknown  function  Y (x) , a circumstance  to 
be  exploited  later  on  in  the  paper. 
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3 • Normalized  System 

The  equations  describing  the  spatial  system  can  be  re- 
written in  a simpler  form  if  the  following  dimensionless 
quantities  are  introduced: 

t=x/L,  m=ML/M^,  y=Y(x)/Y(L),  (18) 

e = a)L/(p/E),  I = M*/M  . (19) 

o 

With  these  definitions,  the  mass  integral  (1)  becomes 


1=  \ mdt  (20) 

J 0 

and  Eqs.  (14)- (16)  become^® 

(d/dt) (my)  + B^my = 0 , (21) 

y(0)  = 0 , (22) 

m(l)y(l)  = 6^  , y(l)  = 1 . (23) 


The  quantity  3,  which  is  proportional  to  the  natural  frequen- 
cy co,  is  called  the  frequency  parameter.  The  additional  bound- 
ary  condition  (23-2)  is  due  to  the  fact  that  the  displacement 
function  Y(x)  has  been  normalized  in  terms  of  Y(L),  the  value 
assumed  by  the  displacement  function  at  the  free  end  of  the 
beam.  This  is  possible,  owing  to  the  fact  that  the  spatial 
system  (14) -(16)  is  homogeneous  in  the  displacement  function 
Y(x). 

1 
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Note  that,  for  problems  where  the  function  m(t)  is  pres- 
cribed a priori  (for  instance,  the  study  of  the  axial  vibration 
of  a constant  section  beam),  the  system  (21) -(23)  is  linear. 

On  the  other  hand,  for  problems  where  the  function  m(t)  is  to 
be  determined  (for  instance,  the  minimum  mass  problem  described 
in  Section  5),  the  system  (21)-(23)  is  nonlinear. 


^^Equation 

placement 


(23-2)  is  a normalization  condition  for  the  dis- 
function y ( t)  . 


■aassT' 


IV 
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4 . Constant-Section  Beam 

In  this  section,  we  consider  the  solution  of  the  system 
(21)- (23)  for  the  particular  case  of  a constant-section  struc- 
ture. For 

m = const  , (24) 

the  mass  integral  (20)  reduces  to 

I = m (25) 

and  Eqs.  (21) -(23)  become 

y + 3^y  = 0 , (26) 

y (0)  = 0 , (27) 

my  (1)  =6^,  y (1)  = 1 . (28) 

The  solution  of  (26),  (27),  and  (28-2),  valid  for  any  value  of 
the  frequency  parameter  S,  is  given  by 

y = sin(6t)/sin6  . (29) 

The  mass  per  unit  length  m corresponding  to  the  frequency 
parameter  6 can  be  computed  from  (28-1) ; 

m = BtanB  . (30) 

Therefore,  the  mass  integral  is  given  by 
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I = 8tanB  . (31) 

In  order  to  understand  the  significance  of  (31) , recall 
that  I is  the  ratio  of  the  beam  mass  to  the  tip  mass  , 
and  consider  two  limiting  cases:  (i)  negligible  mass  ratio  and 
(ii)  infinite  mass  ratio. 

If  1 = 0,  then  the  solution  of  (31)  is 

3 = nTT  , n = 0,1,2, (32) 

On  the  other  hand,  if  I = «>  , then  the  solution  of  (31)  is 

3 = (2n  + 1)tt/2  , n=0,l,2, (33) 

Since  the  first  natural  frequency  or  fundamental  frequency 
corresponds  to  n=  0,  we  conclude  that,  for  any  mass  ratio  in 
the  range 

0 < I < “ , (34) 

the  smallest  value  of  the  frequency  parameter  3 consistent 
with  the  transcendental  equation  (31)  lies  in  the  range 

0 < 3 < Tr/2  . (35) 

We  shall  recall  this  result  in  the  optimization  problem 
formulated  in  Section  5.  Specifically,  we  shall  pose  the  fol- 
lowing problem.  For  a given  value  of  the  frequency  parameter 


16 


AAR-138 


6 in  the  range  (35),  is  there  a beam  yielding  a smaller  value 
of  the  mass  integral  (20)  than  the  value  (31)  associated  with 
the  constant-section  beam?  In  particular,  is  there  a variable- 
section  design  yielding  the  smallest  value  of  the  mass  integral 

(20)  among  all  the  designs  satisfying  the  feasibility  equations 

(21) -(23)? 


1 

■i  ^ 

u 
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5.  Optimization  Problem 

In  the  previous  section,  the  solution  valid  for  a 
constant-section  beam  was  given.  Here,  we  consider  a variable- 
section  beam  and  formulate  the  following  optimization  problem: 
For  a given  value  of  the  frequency  parameter  8,  find  the  mass 
distribution  m(t)  and  the  displacement  distribution  y(t)  such 
that  the  mass  integral  is  minimized,  while  the  feasibility 
equations  are  satisfied.  This  is  a problem  of  the  calculus  of 
variations  or  a problem  of  optimal  control,  depending  on  the 
formulation  being  employed  (Refs.  9-10) . 

Formulation  (FI).  In  this  formulation,  the  optimiza- 
tion problem  has  the  following  format: 

1=  ( mdt  , (36) 

Jo 

my  + my  + 3^niy  = 0 , (37) 

y(0)  = 0 , (38) 

yd)  = 1 , m(l)y(l)  = 3^  . (39) 

Clearly,  the  unknowns  are  the  functions  m(t)  and  y(t). 

Formulation  (F2) . Let  the  following  auxiliary  varia- 
bles be  introduced: 

u=y,  w=y.  (40) 

n 


The  dot  denotes  derivative  with  respect  to  the  independent 
variable  t . 
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With  the  aid  of  these  variables,  the  problem  (36) -(39) 
rewritten  in  optimal  control  format  as  follows: 


y = u, 


m = - (m/u)  (w  + 6 


u = w 


/ 


y(0)  = 0 , 

yd)  = 1 , m(l)u(l)  = 6^  . 


can  be 


(41) 

(42) 

(43) 

(44) 


In  this  formulation,  the  unknowns  are  the  state  variables  y(t), 
m(t),  u(t)  and  the  control  variable  w(t)  . Note  that  the  control 
appears  linearly  in  the  differential  equations.  Hence,  the 
Hamiltonian  is  linear  in  the  control  w(t),  and  the  problem  is 
singular.  Also  note  that  the  system  (41) -(44)  is  autonomous, 
since  the  independent  variable  t does  not  appear  explicitly  on 
the  right-hand  side  of  Eqs.  (41)-(42).  As  a consequence,  the 
Hamiltonian  is  constant  along  the  interval  of  integration. 


I 
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6 . Alternative  Formulations  of  the  Optimization  Problem 
In  this  section,  we  present  some  alternative  formula- 
tions of  the  optimization  problem,  having  important  computa- 
tional implications. 

Formulation  (F3) . Let  the  following  auxiliary  variable 
be  introduced: 

z = mu  . (45) 

Let  the  variable  m(t)  of  Formulation  (F2)  be  eliminated 

and  replaced  by  z(t).  With  this  understanding,  the  problem 
(41)- (44)  can  be  reformulated  as  follows: 


I = \ (z/u)dt  , 

(46) 

•'0 

u , z = -B^yz/u  , u = w , 

(47) 

y(0)  = 0 , 

(48) 

yd)  = 1 , z(l)  = . 

(49) 

In  this  formulation,  the  unknowns  are  the  state  variables  y(t), 
z(t),  u(t)  and  the  control  variable  w(t) . After  these  functions 
are  determined,  the  function  m(t)  can  be  computed  a posteriori 
with  the  relation 

m = z/u  . (50) 


I 
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A characteristic  of  this  formulation  is  that  all  the 
boundary  conditions  are  linear.  The  Hamiltonian  is  linear  in 
the  control  w(t),  and  hence  the  problem  is  singular.  Once 
more,  the  system  is  autonomous,  and  the  Hamiltonian  is  cons- 
tant along  the  interval  of  integration. 

Formulation  (F4).  Inspection  of  (46) -(49)  shows  that 
the  control  w(t)  appears  only  in  Eq.  (47-3).  Also,  both 
u(0)  and  u(l)  are  free.  Hence,  one  surmises  that  the  differential 
constraint  (47-3)  might  be  redundant.  This  can  be  shown  as  follows. 
Assume  that  the  problem  (46) -(49)  is  solved  bypassing  Eq. 

(47-3);  with  the  solution  y(t),  z(t),  u(t)  known,  one  can  al- 
ways determine  a posteriori  the  function  w(t)  in  such  a way 
that  the  feasibility  equation  (47-3)  is  satisfied.  This  simple 
but  important  observation  leads  to  a new  formulation  of  the 
optimal  control  problem: 


I = (z/u)dt  , (51) 

Jo 

y = u , z = -6^yz/u  , (52) 

y(0)  = 0 , (53) 

yd)  = 1 , z(l)  = . (54) 


i 


In  this  formulation,  the  unknowns  are  the  state  variables  y(t), 
z(t)  and  the  control  variable  u(t).  After  these  functions  are 
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determined,  the  functions  m(t)  and  w(t)  are  computed  a post- 
eriori with  the  relations 


m = z/u  , 


A characteristic  of  this  formulation  is  that  all  the 
boundary  conditions  are  linear.  In  addition,  the  state  varia- 
ble u(t)  of  Formulation  (F3)  becomes  the  control  variable 
u(t)  of  Formulation  (F4).  The  Hamiltonian  is  nonlinear 
in  the  control  u(t),  and  hence  the  problem  is  nonsingular. 

Once  more,  the  system  is  autonomous,  and  the  Hamiltonian  is 
constant  along  the  interval  of  integration. 

Formulation  (F5) . This  formulation  is  a slight  modifi- 
cation of  the  previous  one  and  is  obtained  by  eliminating  the 
control  variable  u(t)  and  replacing  it  with  the  new  control 
variable  m(t).  On  account  of  (55-1),  Eqs.  (51) -(54)  become 


I = l mdt  , 


y = z/m  , 


z = -6  ym  , 


y(0)  = 0 , 


yd)  = 1 , 


z(i)  = e' 


In  this  formulation,  the  unknowns  are  the  state  variables  y(t), 
z(t)  and  the  control  variable  m(t).  After  these  functions  are 
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determined,  the  functions  u(t)  and  w(t)  are  computed  a post- 
eriori with  the  relations 

u = z/m  , w = u . (60) 

A characteristic  of  this  formulation  is  that  all  the 
boundary  conditions  are  linear.  In  addition,  the  Hamiltonian 
is  nonlinear  in  the  control  m(t) , and  hence  the  problem  is  non- 
singular. Once  more,  the  system  is  autonomous,  and  the  Hamil- 
tonian is  constant  along  the  interval  of  integration.  Note 
that  this  formulation  can  be  obtained  directly  from  (20) -(23) 
after  consideration  of  the  relation 

z = my  . (61) 

Remark.  The  analytical  justification  for  the  passage 
from  Formulation  (F3)  to  Formulation  (F4)  is  as 
follows.  With  reference  to  Formulation  (F3),  let  Xj^(t), 

A2(t),  X2(t)  denote  variable  Lagrange  multipliers  associated 
with  the  differential  constraints  (47)  . Let 

H = z/u  - X j^u  + w (62) 


denote  the  Hamiltonian  of  the  system  (46)-(49).  The  condition 
optimizing  the  control  w(t)  of  Formulation  (F3)  takes  the 
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form 


= 0 ' 0 < t < 1 , (63) 

which  implies  that 

A3 (t)  = 0 , 0 < t < 1 . (64) 

Note  that  both  u{0)  and  u(l)  are  free;  hence,  the  transversality 
condition  requires  that 

A3{0)  = A3(1)  = 0 , (65) 

and  these  relations  are  consistent  with  (64).  From  (64)-(65), 
it  is  natural  to  infer  that  the  differential  constraint  (47-3) 
is  redundant,  a circumstance  which  has  important  analytical 
and  computational  implications. 
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7 . Description  of  the  Algorithms 

With  the  sequential  ordinary  gradient-restoration  al- 
gorithm (SOGRA,  Refs.  3-4),  numerical  solutions  can  be  obtained 
using  any  of  the  optimal  control  formulations  of  Sections  5-6. 
This  is  because  SOGRA  can  be  applied  regardless  of  whether  the 
problem  is  singular  or  nonsingular.  On  the  other  hand,  with 
the  modified  quasilinearization  algorithm  (MQA,  Refs.  5-6), 
numerical  solutions  can  be  obtained  using  Formulations  (F4) 
and  (F5) , but  not  Formulations  (F2)  and  (F3) . This 
is  because  the  optimal  control  problem  is  singular  with  the 
latter  formulations,  while  it  is  nonsingular  with  the  former. 

For  these  reasons,  the  experiments  reported  here  refer  to 
Formulations  (F4)  and  (F5) . 

Note  that  both  SOGRA  and  MQA  require  the  state  vector  to  be 
Icnown  at  the  initial  point,  so  that  the  trajectories  in  the 
state  space  can  be  anchored  at  the  initial  point.  With  refer- 
ence to  the  fourth  and  fifth  formulations  (state  vector  given 
at  the  final  point) , this  situation  can  be  achieved  by  re- 
placing the  independent  variable  t with  its  complement 

T = 1-t  , (66) 

that  is,  by  reversing  tlie  sense  of  integration  and  inter- 
changing the  final  point  with  the  initial  point.  Even  though 
this  replacement  of  independent  variable  is  necessary 
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computationally,  it  is  not  formally  executed  here,  for  two 
reasons;  for  the  sake  of  brevity  and  in  order  to  avoid  compli- 
cating the  description  of  the  problems  and  the  interpretation 
of  the  results.  Therefore,  in  the  following  sections,  we 
shall  retain  the  independent  variable  t,  that  is,  we  shall 
imagine  the  trajectories  in  the  state  space  to  be  anchored  at 
the  final  point  at  every  iteration. 

Formulation  (F4).  With  the  formulation  represented  by 

12 

Eqs.  (51) -(54),  the  Hamiltonian  is  given  by 

H = z/u  - Xj^u  + A28^yz/u  , (67) 


where  Xj^(t)  and  X2(t)  denote  variable  Lagrange  multipliers 

associated  with  the  differential  constraints  (52) . Hence,  the 

12 

optimality  conditions  take  the  form 


X^  + (z/u^)  (1  + X26^y)  = 0 , 

(68) 

= X2  = (1/u)  (1  + X26^y)  , 

(69) 

Xj^(O)  + u = 0 , X2(0)  = 0 . 

(70) 

12 

Equation  (68)  optimizes  the  control  distribution 
(70-1),  the  symbol  y denotes  a constant  Lagrange 
associated  with  the  initial  condition  (53) . 

u (t) . In  Eq 
multiplier 
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Formulation  (F5) . With  the  formulation  represented 
by  Eqs.  (56) -(59),  the  Hamiltonian  is  given  by 

2 

H = m - A ^z/m + A 23  ym  , (71) 

where  Aj^(t)  and  A2(t)  denote  variable  Lagrange  multipliers 

associated  with  the  differential  constraints  (57) . Hence,  the 

13 

optimality  conditions  take  the  form 

Aj^z/m^  + 1 + A23^y  = 0 , (72) 

• 2 

^1  ~ ^2^  m , A2=-Aj^/m  , (73) 

A^ (0)  + M = 0 , A2 (0)  = 0 . (74) 


Auxiliary  Functionals.  In  addition  to  the  functional  I, 
several  auxiliary  functionals  are  of  interest  in  the  implemen- 
tation of  SOGRA  and  MQA,  These  auxiliary  functionals  are 
denoted  by  the  symbols  J,  P,  Q,  R and  have  the  following  mean- 
ing: J is  the  augmented  functional,  that  is,  the  functional  I 


Equation  (72)  optimizes  the  control  distribution  ro(t).  In 
Eq.  (74-1),  the  symbol  y denotes  a constant  Lagrange  multi- 
plier associated  with  the  initial  condition  (58). 
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augmented  linearly  by  the  constraints  through  Lagrange  multi- 
pliers; P is  the  constraint  error,  that  is,  the  norm  squared 
of  the  error  in  the  feasibility  equations;  Q is  the  optimality 
condition  error,  that  is,  the  norm  squared  of  the  error  in  the 
optimality  conditions;  and  R=P+Q  is  the  total  error  in  the 
system  composed  of  the  feasibility  equations  and  the  optimali- 
ty conditions.  In  the  definitions  of  J,  P,  Q,  R,  it  is  tacitly 
assumed  that  the  final  conditions  (54)  or  (59)  are  satisfied  at 
every  iteration  of  SOGRA  and  MQA. 

Convergence  Conditions.  The  auxiliary  functionals  P,  Q, 

R are  defined  in  such  a way  that 

P = Q=R=0  (75) 

for  the  optimal  solution.  For  an  approximation  to  the  optimal 
solution,  Eqs.  (75)  are  replaced  by  the  inequalities 

P<Cj^,  Q < , (76) 

or 

R<e3,  (77) 

where  £2'  ^3  small,  preselected  numbers.  Inequalities 

(76)  are  of  interest  for  SOGRA,  and  Ineq.  (77)  is  of  interest 


for  MQA. 
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Sequential  Ordinary  Gradient-Restoration  Algorithm.  The 
sequential  ordinary  gradient-restoration  algorithm  (SOGRA, 

Refs.  3-4)  is  an  iterative  technique  which  includes  a sequence 
of  cycles  having  the  following  properties:  (i)  the  functions 
available  both  at  the  beginning  and  at  the  end  of  each  cycle 
are  feasible;  that  is,  they  are  consistent  with  the  feasibi- 
lity equations  within  the  preselected  accuracy  (76-1);  and  (ii) 
the  functions  produced  at  the  end  of  each  cycle  are  charac- 
terized by  a value  of  the  functional  I which  is  smaller  than 
that  associated  with  the  functions  available  at  the  beginning 
of  the  cycle. 

To  achieve  the  above  properties,  each  cycle  is  made  of 
two  phases,  a gradient  phase  and  a restoration  phase.  These 
phases  are  now  described. 

The  gradient  phase  is  started  only  when  the  constraint 
error  P satisfies  Ineq.  (76-1) . It  involves  a single  iter- 
ation, which  is  designed  to  decrease  the  value  of  the  func- 
tional I or  the  augmented  functional  J,  while  satisfying  the 
constraints  to  first  order.  During  this  iteration,  the  first 
variation  of  the  functional  I is  minimized,  subject  to  the 
linearized  constraints  and  a quadratic  constraint  on  the  vari- 
ations of  the  control. 

The  restoration  phase  is  started  only  when  the  constraint 
error  P violates  Ineq.  (76-1) . The  restoration  phase  involves 
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one  or  more  iterations.  In  each  restorative  iteration,  the 
objective  is  to  reduce  the  constraint  error  P,  while  the  con- 
straints are  satisfied  to  first  order  and  the  norm  of  the 
variations  of  the  control  is  minimized.  The  restoration  phase 
is  terminated  whenever  Ineq.  (76-1)  is  satisfied. 

In  summary,  the  main  properties  of  the  sequential  ordi- 
nary gradient-restoration  algorithm  can  be  written  as  follows: 

J < J , P < f (78) 

P < P , P > , (79) 

i < I , P < , P 1 , (80) 

where  (78)  hold  for  a gradient  iteration,  (79)  hold  for  a 
restorative  iteration,  and  (80)  hold  for  a complete  gradient- 
restoration  cycle.  In  the  above  relations,  I,  J,  P denote 
quantities  evaluated  at  the  beginning  of  an  iteration  or  a 
cycle,  while  I,  J,  P denote  quantities  evaluated  at  the  end 
of  an  iteration  or  a cycle. 

At  each  iteration  of  the  gradient  phase  or  the  restora- 
tion phase,  a linear,  two-point  boundary-value  problem  must  be 
solved.  If  the  method  of  particular  solutions  is  employed 
(Refs.  11-12),  one  must  execute  q+1  independent  sweeps  of  the 
linearized  system  governing  the  variations  associated  with  the 

f 

1 

I 
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gradient  phase  or  the  restoration  phase,  where  q is  the  number 
of  initial  conditions.  Since  q=  1 for  the  present  problem,  each 
gradient  iteration  or  restorative  iteration  requires  2 sweeps 
of  the  linearized  system.  For  the  details,  see  Refs.  3-4  and 
11-12. 

Modified  Quasilinearization  Algorithm.  The  modified 
quasilinearization  algorithm  (MQA,  Refs.  5-6)  is  an  iterative 
technique  which  includes  a sequence  of  iterations  having  the 
following  property;  the  functions  produced  at  the  end  of  each 
iteration  are  characterized  by  a value  of  the  functional  R 
which  is  smaller  than  that  associated  with  the  functions  avail- 
able at  the  beginning  of  the  iteration.  Therefore,  the  main 
property  of  the  modified  quasilinearization  algorithm  can  be 
written  as  follows: 

R<R,  (81) 

with  the  implication  that 

P < P and/or  Q < Q . (82) 

At  each  iteration,  both  the  feasibility  equations  and 
the  optimality  conditions  are  satisfied  to  first  order.  To 
achieve  this,  a linear,  two-point  boundary-value  problem  must 
be  solved.  If  the  method  of  particular  solutions  is  employed 
(Refs.  11-12),  one  must  execute  n + 1 independent  sweeps  of  the 
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linearized  system  governing  the  variations  associated  with 
MQA,  where  n is  the  number  of  state  variables.  Since  n = 2 for 
the  present  problem,  each  iteration  requires  3 sweeps  of  the 
linearized  system.  For  the  details,  see  Refs.  5-6  and  11-12. 


1 

I 


32 


AAR-138 


8 . Experimental  Conditions 

The  problem  of  determining  the  mass  distribution  that 
minimizes  the  total  mass  for  a given  value  of  the  frequency 
parameter  3 was  solved  using  both  Formulation  (F4)  and 
Formulation  (F5) . The  following  values  were  given  to  the 
frequency  parameter: 

B=Tr/6,  3 = tt/5  , g = Tr/4  , 3 = tt/3  , 3 = tt/2  . (83) 

Computations  were  performed  on  the  IBM  370/155  computer 
of  Rice  University,  Houston,  Texas,  using  both  the  sequential 
ordinary  gradient-restoration  algorithm  (SOGRA)  and  the  modi- 
fied quasilinearization  algorithm  (MQA) . Both  algorithms  were 
programmed  in  FORTRAN  IV,  and  the  numerical  results  were  ob- 
tained in  double-precision  arithmetic. 

The  interval  of  integration  was  divided  into  100  steps 
for  SOGRA  and  50  steps  for  MQA.  The  differential  eauations  were 
integrated  using  Hamming's  modified  predictor-corrector  method, 
with  a special  Runge-Kutta  starting  procedure.  The  definite  inte- 
grals I,  J,  P,  Q were  computed  using  a modified  Simpson's 
rule. 

Nominal  Functions.  For  each  formulation,  two  sets  of 
nominal  functions  were  employed,  more  specifically: 


(Nl)  y (t)  = t , 


z(t)  = 3^, 


u(t)  = 1 , 


(84) 
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(N2) 

y (t)  = t^  , 

z (t)  = 

u (t)  = 1 , 

(85) 

and 

(N3) 

y(t)  = t , 

z ( t)  = 

6^  , 

m ( t)  =1  , 

(86) 

(N4) 

Y(t)  = t^. 

z (t)  = 

m (t)  = 1 . 

(87) 

The 

nominal  functions 

(84)  and 

(85)  were  employed 

in  connec- 

tion 

with 

Formulation 

(F4)  , 

and 

the  nominal 

functions 

(86)  and  (87)  were  employed  in  connection  with  Formulation 
(F5) . Note  that  these  nominal  functions  satisfy  the  bound- 

ary conditions,  but  violate  the  differential  constraints. 

Stopping  Conditions.  The  sequential  ordinary  gradient- 
restoration  algorithm  was  programmed  to  stop  whenever  a 
solution  consistent  with  the  following  inequalities  was 
obtained : 

P < E-08  , Q < E-04  . (88) 

The  modified  quasilinearization  algorithm  was  programmed  to 
stop  whenever  the  following  inequality  was  satisfied: 

R < E-  10  , (89) 


12 

The  symbol  E + ab  stands  for  10 


I 

I 
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with  the  implication  that 

P < E-10  , Q < E-10.  (90) 

Remark  3.1.  For  the  sake  of  brevity,  the  computational 
details  pertaining  to  SOGRA  and  MQA  (for  instance,  solution  of 
the  linear,  two-point  boundary  value  problem,  computation  of 
the  stepsize,  safeguards,  and  nonconvergence  conditions)  are 
omitted,  and  the  reader  is  referred  to  Refs.  3-4  and  Refs. 

5-6. 

Remark  8.2.  Ordinarily,  algorithms  of  the  second-order 
type  require  the  specification  of  nominal  functions  not  only 
for  the  state  and  the  control,  but  also  for  the  multipliers. 

This  is  not  the  case  with  MQA.  Once  the  nominal  state  and  the 
nominal  control  are  specified,  the  nominal  multipliers  X^(t), 
A2{t),  y are  determined  automatically  through  a subroutine 
which  supplies  the  solution  of  the  following  auxiliary  minimi- 
zation problem;  minimize  the  error  in  the  optimality  conditions 
Q with  respect  to  multipliers  Aj^(t),  ^2^^),  y-  This  is  an 
important  feature  of  MQA,  certainly  very  useful  to  engineers, 
since  it  circumvents  the  need  for  an  arbitrary  guess  of  the 
multiplier  functions. 
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9 . Numerical  Results 

Under  the  experimental  conditions  outlined  in  the  pre- 
vious section,  40  computer  runs  v?ere  executed:  the  minimization 
problem  was  solved  for  five  different  values  of  the  frequency 
parameter,  using  two  algorithms  (SOGRA  and  MQA) , two  formula- 
tions [(F4)  and  (F5)],  and  two  sets  of  nominal  functions  for 
each  formulation  [(Nl)  and  (N2)  for  Formulation  {F4),  and  (N3) 
and  (N4)  for  Formulation  (F5)].  Convergence  to  the  desired 
stopping  condition  was  achieved  in  37  runs.  The  nonconvergence 
occurring  in  the  remaining  3 runs  was  due  to  poor  choice  of  the 
nominal  functions.  Incidentally,  all  the  nonconvergence  cases 
occurred  for  3 = fr/2 , which  is  the  upper  limit  to  the  range  of 
values  of  the  frequency  parameter  investigated.  The  numerical 
results  are  presented  in  Tables  1-24. 

Tables  1-4  present  summary  results  pertaining  to  SOGRA 
at  convergence.  For  each  value  of  the  frequency  parameter  3, 
the  tables  show  the  number  of  iterations  for  convergence  N, 
the  computed  value  of  the  functional  I,  the  exact  value  of  the 
functional  , the  number  of  correct  significant  digits  M 
(determined  by  comparing  I with  I^) , and  the  error  in  the 
optimality  conditions  Q.  Clearly,  as  far  as  the  minimum  mass 
is  concerned,  the  solutions  obtained  are  precise  to  at  least 
4 significant  digits;  in  some  cases,  they  are  precise  to  6 
significant  digits. 
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Tables  5-8  present  summary  results  pertaining  to  MQA 
at  convergence.  For  each  value  of  the  frequency  parameter  6, 
the  tables  show  the  number  of  iterations  for  convergence  N, 
the  computed  value  of  the  functional  I,  the  exact  value  of 
the  functional  the  number  of  correct  significant  digits 

M (determined  by  comparing  I with  I^) , and  the  total  error 
in  the  system  R.  Clearly,  as  far  as  the  minimum  mass  is  con- 
cerned, the  solutions  obtained  are  precise  to  at  least  6 
significant  digits;  in  many  cases,  they  are  precise  to  7 sig- 
nificant digits. 

Due  to  the  large  number  of  runs  performed,  it  is 
impractical  to  present  the  convergence  history  for  every  value 
of  the  frequency  parameter  B.  However,  for  a particular  value 
of  the  frequency  parameter  (namely,  B = Tr/4),  the  convergence 
history  of  SOGRA  is  given  in  Tables  9-12  and  the  convergence 
history  of  MQA  is  given  in  Tables  13-16.  The  descent  property 
in  I,  characteristic  of  SOGRA,  is  apparent  from  Tables  9-12. 
Analogously,  the  descent  property  in  R,  characterisitc  of 
MQA,  is  apparent  from  Tables  13-16. 

Due  to  the  large  number  of  runs  performed,  it  is  imprac- 
tical to  present  the  converged  solutions  y(t),  z(t),  m(t),  u(t), 
w(t)  for  every  case.  Therefore,  Tables  17-21  present  only  one 
set  of  converged  solutions;  these  are  the  solutions  corres- 
ponding to  the  lowest  value  of  R achieved  computationally  for 
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any  given  value  of  the  frequency  parameter  3 . Even  though 
the  tolerance  level  set  for  MQA  is  R<E-10,  values  of  the 
total  error  much  lower  that  E-10  were  achieved  during  compu- 
tation, more  precisely,  R < E-16  for  3 = fr/S  and  R < E-18  for 
the  remaining  values  of  3. 

Table  22  is  obtained  from  Tables  17-21  by  normalizing 
the  mass  distribution  m(t)  with  respect  to  the  mass  per  unit 
length  at  t=0.  Therefore,  Table  22  shows  the  ratio  m(t)/m(0) 
for  several  values  of  the  frequency  parameter  3.  Clearly, 
the  mass  distribution  of  the  optimal  beam  approaches  that  of 
the  constant-section  beam  for  3 0 and  deviates  considerably 

from  it  for  3^fT/2.  In  particular,  the  ratio  m{l)/m{0)  of 
tip  mass  per  unit  length  to  root  mass  per  unit  length  approa- 
ches 1 for  3 -»■  0 and  1/6  for  3->-tt/2. 

Finally,  Tables  23-24  show  the  mass  integral  for  the 
optimal  beam  I^,  the  mass  integral  for  the  constant-section 
beam  I , as  well  as  the  quantities  I /I  and  (I  - I )/I  for 
several  values  of  the  frequency  parameter  3.  For  3 -*■  0 , the 
optimal  beam  exhibits  no  weight  advantage  over  the  constant- 
section  beam.  For  3=  0-5,  the  weight  advantage  is  only  0.6%. 
For  3=1.0,  the  weight  advantage  is  11.3%.  For  3 = 1.4,  the 
weight  advantage  is  55.3%.  Finally,  for  the  weight 

advantage  approaches  100%.  The  explanation  for  this  limiting 
result  is  simple;  for  the  mass  of  the  optimal  beam  is 

finite,  while  the  mass  of  the  constant-section  beam  is  infi- 
nitely large. 
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Table  1.  Summary  results,  SOGRA, 

Formulation  (F4),  Nominal  (Nl). 


8 

N 

I 

I 

e 

M 

Q 

7T/6 

3 

0.3001435E+00 

0.3001434E+00 

6 

0.85E-07 

tt/S 

3 

0.4495492E+00 

0.4495488E+00 

5 

0.85E-06 

jr/4 

3 

0.7545937E+00 

0.7545892E+00 

4 

0.15E-04 

tt/3 

4 

0.1560903E+01 

0.1560918E+01 

5 

0.28E-04 

7T/2 

10 

0.5295951E+01 

0. 5295977E+01 

5 

0.88E-03  (*) 

(*)  Algorithm  unable  to  reach  desired  stopping  condition; 

loss  of  descent  property  on  I due  to  numerical  inaccuracy . 


Table  2.  Summary  results,  SOGRA, 

Formulation  (F4),  Nominal  (N2). 


3 

N 

I 

I 

e 

M 

0 

■it/6 

4 

0.3001435E+00 

0.3001434E+00 

6 

0.85E-07 

■n/5 

4 

0.4495492E-i-00 

0.4495488E+00 

5 

0.85E-06 

Tr/4 

4 

0.7545937E+00 

0.7545892E+00 

4 

0.15E-04 

tt/3 

5 

0.1560903E4-01 

0.1560918E+01 

5 

0.28E-04 

ti/2 

11 

0.5295951E+01 

0.5295977E+01 

5 

0.88E-03(*) 

{*)  Algorithm  unable  to  reach  desired  stopping  condition; 

loss  of  descent  property  on  I due  to  numerical  inaccuracy. 


y 
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Table  3.  Summary  results,  SOGRA, 

Formulation  (F5),  Nominal  (N3). 


6 

N 

T 

^e 

M 

Q 

7t/6 

19 

0.3001453E+00 

0.3001434E+00 

5 

0. 32E-04 

it/5 

14 

0.4495573E+00 

0.4495488E+00 

4 

0.84E-04 

7t/4 

10 

0.7545631E+00 

0.7545892E+00 

4 

0.36E-04 

tt/3 

12 

0.1560911E+01 

0.1560918E+01 

6 

0. 51E-04 

tt/2 

(*) 

(*) 

0.5295977E+01 

(*) 

(*) 

(*)  Nonconvergence,  failure  of  initial  restoration  phase; 
algorithm  unable  to  find  a feasible  solution  in  less 
than  10  restorative  iterations. 


Table  4.  Summary  results,  SOGRA, 

Formulation  (F5),  Nominal  (N4). 


B 

N 

I 

I 

e 

M 

Q 

7I/6 

21 

0.3001449E+00 

0.3001434E+00 

5 

0. 25E-04 

•it/5 

14 

0.4495525E+00 

0.4495488E+00 

4 

0.38E-04 

it/4 

10 

0.7545679E+00 

U.7545892E+00 

4 

0.26E-04 

Tr/3 

15 

0.1560913E+01 

0.1560918E+01 

6 

0. 54E-04 

tt/2 

(*) 

(*) 

0. 5295977E+01 

(*) 

(*) 

{*)  Nonconvergence,  failure  of  initial  restoration  phase; 
algorithm  unable  to  find  a feasible  solution  in  less 
than  10  restorative  iterations. 
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Table  5-  Summary  results,  MQA, 

Formulation  (F4),  Nominal  (Nl). 


N 

I 

I 

e 

M 

R 

u/ 6 

3 

0.3001434E+00 

0.3001434E+00 

7 

0.20E-18 

it/5 

3 

0.4495488E+00 

0.4495488E+00 

7 

0.74E-16 

tt/4 

3 

0.7545892E+00 

0.7545892E+00 

7 

0.12E-12 

Tr/3 

4 

0.1560918E+01 

0.1560918E+01 

7 

0.32E-18 

Tr/2 

(*) 

(*) 

0.5295977E+01 

(*) 

(*) 

(*)  Nonconvergence,  stepsize 

bisection  limit 

reached . 

Table  6.  Summary  results,  MQA, 

Formulation  (F4),  Nominal 

(N2)  . 

6 

N 

I 

I 

e 

M 

R 

Tr/6 

3 

0.3001434E+00 

0.3001434E+00 

7 

0.17E-18 

tt/5 

3 

0.4495488E+00 

0.4495488E+00 

7 

0.70E-16 

7r/4 

3 

0.7545892E+00 

0.7545892E+00 

7 

0.97E-13 

■Tr/3 

4 

0.1560918E+01 

0.1560918E+01 

7 

O.llE-18 

7r/2 

6 

0.5295977E+01 

0.5295977E+01 

7 

0.15E-15 
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Table  7.  Summary  results,  MQA, 

Formulation  (F5),  Nominal  (N3). 


3 

N 

I 

I 

e 

M 

R 

tt/6 

9 

0.3001433E+00 

0.3001434E+00 

6 

0.79E-10 

tt/5 

6 

0.4495488E+00 

0.4495488E+00 

7 

0.15E-14 

77/4 

5 

0.7545892E+00 

0.7545892E+00 

7 

0.41E-18 

77/3 

4 

0.1560917E+01 

0.1560918E+01 

6 

0. 20E-10 

77/2 

7 

0. 5295976E+01 

0.5295977E+01 

6 

0.38E-17 

Table  8.  Summary  results,  MQA, 

Formulation  (F5),  Nominal 

(N4)  . 

6 

N 

I 

I 

e 

M 

R 

77/6 

9 

0.3001433E+00 

0.3001434E+00 

6 

0.59E-10 

77/5 

7 

0.4495487E+00 

0.4495488E+00 

6 

0. 55E-10 

7t/4 

4 

0.7545890E+00 

0.7545892E+00 

6 

0.88E-11 

77/3 

4 

0.1560917E+01 

0.1560918E+01 

6 

0.83E-10 

77/2 

7 

0.5295976E+01 

0. 5295977E+01 

6 

0.78E-18 

42 


AAR-138 


Table 

9.  Convergence  history,  3 = tt/4  , SOGRA, 
Formulation  (F4),  Nominal  (Nl). 

N 

Phase 

P 

Q 

I 

0 

0.48E-01 

1 

REST 

O.lOE-33 

0.21E-01 

0.76081 

2 

GRAD 

0.13E-04 

3 

REST 

0.82E-33 

0.15E-04 

0.75459 

Table 

10.  Convergence  history,  3 = Tr/4 , SOGRA, 

Formulation  (F4),  Nominal  (N2) . 

N 

Phase 

P 

Q I 

0 

0.36E+00 

1 

REST 

0.29E-03 

2 

REST 

0.90E-33 

0.21E-01  0.76081 

3 

GRAD 

0.13E-04 

4 


REST 


0. 27E-32 


0.15E-04 


0.75459 
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Table  11.  Convervenge  history,  6 = Tr/4 , SOGRA, 
Formulation  (F5) , Nominal  (N3). 


N 

Phase 

P 

Q 

I 

0 

0,27E+00 

1 

REST 

0.13E+00 

2 

REST 

O.llE-01 

3 

REST 

0.58E-04 

4 

REST 

0.37E-08 

0.31E+00 

0.83009 

5 

GRAD 

0.16E-01 

6 

REST 

0.80E-04 

7 

REST 

0.41E-08 

0.46E-02 

0.75504 

8 

GRAD 

0.25E-05 

9 

REST 

0. 26E-12 

0.24E-03 

0.75462 

10 

GRAD 

0.95E-08 

0.36E-04 

0.75456 

Table  12.  Convergence 
Formulation 

history,  3 = ■'r/4  , SOGRA, 
(F5) , Nominal  (N4). 

N Phase  P 

Q 

I 

0 

0. 55E+00 

1 

REST 

0.81E-01 

2 

REST 

0.38E-02 

3 

REST 

0.92E-05 

4 

REST 

O.lOE-09 

0. 28E+00 

0.82212 

5 

GRAD 

0.13E-01 

6 

REST 

0.48E-04 

7 

REST 

0.15E-08 

0.46E-02 

0.75503 

8 

GRAD 

0.27E-05 

9 

REST 

0.20E-12 

0.19E-03 

0.75461 

10 

GRAD 

0. 56E-08 

0.26E-04 

0.75456 
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Table  13.  Convergence  history,  B = tt/4  , MQA, 
Formulation  (F4) , Nominal  (Nl) . 


N 

P 

Q 

R 

I 

0 

0.48E-01 

0.32E-01 

0.80E-01 

0.61685 

1 

0.84E-05 

0.23E-02 

0.23E-02 

0.75463 

2 

0.48E-08 

0. 96E-06 

0.96E-06 

0.75455 

3 

0.44E-15 

0.12E-12 

0.12E-12 

0.75458 

Table 

14.  Convergence  history,  B = tt/4,  MQA, 
Formulation  (F4),  Nominal  (N2). 

N 

P 

Q 

R 

I 

0 

0.36E+00 

0.27E-01 

0.39E+00 

0.61685 

1 

0.27E-03 

0.30E-02 

0.33E-02 

0.74916 

2 

0.72E-08 

0.75E-06 

0.76E-06 

0.75454 

3 

O.lOE-14 

0.96E-13 

0.97E-13 

0.75458 
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Table  15.  Convergence  history,  3 = ''t/4  , MQA, 
Formulation  (F5) , Nominal  (N3). 


N 

P 

Q 

R 

I 

0 

0. 27E+00 

0.20E+00 

0.47E+00 

1.00000 

1 

0.48E-01 

0. 94E-01 

0.14E+00 

0.75124 

2 

0.30E-03 

0.71E-02 

0.74E-02 

0.74933 

3 

0.69E-06 

0.26E-04 

0.27E-04 

0.75419 

4 

0. 21E-10 

0, 58E-09 

0.61E-09 

0.75458 

5 

0.15E-19 

0.39E-18 

0.41E-18 

0.75458 

Table 

16.  Convergence  history,  3 = Tr/4,  MQA, 
Formulation  (F5) , Nominal  (N4). 

N 

P 

Q 

R 

I 

0 

0. 55E+00 

0.16E+00 

0.72E+00 

1.00000 

1 

0. 98E-01 

0.45E-01 

0.14E+00 

0.77369 

2 

0.21E-03 

0.19E-02 

0.22E-02 

0.74914 

3 

0.14E-06 

0.29E-05 

0.30E-05 

0.75442 

4 

0.33E-12 

0.84E-11 

0.88E-11 

0.75458 

I 

I 

I 
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Table  17. 

Converged 

solution. 

B = 7r/6,  MQA. 

t 

y 

z 

m 

u 

w 

0.0 

0.0000 

0.3126 

0.3270 

0.9557 

0.0000 

0.1 

0.0956 

0.3121 

0.3261 

0.9570 

0.0262 

0.2 

0.1914 

0.3108 

0.3235 

0.9609 

0.0524 

0.3 

0.2878 

0.3087 

0.3191 

0.9675 

0.0789 

0.4 

0.3850 

0.3058 

0.3131 

0.9767 

0.1055 

0.5 

0.4833 

0.3021 

0.3056 

0.9886 

0.1325 

0.6 

0.5829 

0.2977 

0.2968 

1.0032 

0.1598 

0.7 

0.6840 

0.2927 

0.2868 

1.0206 

0.1875 

0.8 

0.7871 

0.2870 

0.2757 

1.0408 

0.2157 

0.9 

0.8923 

0.2808 

0.2639 

1.0638 

0.2446 

1.0 

1.0000 

0.2741 

0.2515 

1.0897 

0.2741 
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Table  18. 

Converged 

solution , 3 = 

Tr/5,  MQA. 

t 

y 

z 

m 

u 

w 

0.0 

0.0000 

0.4753 

0.5072 

0.9371 

0.0000 

0.1 

0.0937 

0.4743 

0.5052 

0.9389 

0.0370 

0.2 

0.1879 

0.4715 

0.4992 

0.9445 

0.0741 

0.3 

0.2828 

0.4669 

0.4896 

0.9538 

0.1116 

0.4 

0.3788 

0.4606 

0.4764 

0.9668 

0.1495 

0.5 

0.4763 

0.4527 

0.4602 

0.9837 

0.1880 

0.6 

0.5756 

0.4434 

0.4414 

1.0044 

0.2272 

0.7 

0.6773 

0.4327 

0.4204 

1.0292 

0.2674 

0.8 

0.7816 

0.4209 

0.3979 

1.0580 

0.3085 

0.9 

0.8890 

0.4082 

0.3742 

1.0909 

0.3509 

1.0 

1.0000 

0.3947 

0.3499 

1.1282 

0.3947 

Table  19. 

48 

Converged 

solution, 

AAR- 

6 = tt/4,  MQA. 

•138 

t 

y 

z 

m 

u 

w 

0.0 

0.0000 

0.8170 

0.9037 

0.9041 

0.0000 

0.1 

0.0905 

0.8145 

0.8981 

0.9069 

0.0558 

0.2 

0.1815 

0.8071 

0.8817 

0.9153 

0.1120 

0.3 

0.2737 

0.7949 

0.8553 

0.9293 

0.1688 

0.4 

0.3676 

0.7783 

0.8200 

0.9491 

0.2267 

0.5 

0.4637 

0.7578 

0.7775 

0.9747 

0.2860 

0,6 

0.5627 

0.7340 

0.7293 

1.0063 

0.3471 

0.7 

0.6652 

0.7074 

0.6774 

1.0442 

0.4103 

0.8 

0.7718 

0.6786 

0.6234 

1.0885 

0.4761 

0.9 

0.8832 

0.6482 

0.5688 

1.1395 

0.5448 

1.0 

1.0000 

0.6168 

0.5150 

1.1976 

0.6168 
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Table  20.  Converged  solution,  B = it/3  , MQA. 


z m u w 


o 

o 

0.0000 

1.7549 

2.0937 

0.8381 

0.0000 

0.1 

0.0839 

1.7453 

2.0709 

0.8427 

0.0920 

0.2 

0.1688 

1.7171 

2.0044 

0.8566 

0.1851 

0.3 

0.2556 

1.6717 

1.8999 

0.8798 

0.2803 

0.4 

0.3451 

1.6114 

1.7654 

0.9127 

0.3785 

0.5 

0.4385 

1.5390 

1.6103 

0.9557 

0.4808 

0.6 

0.5366 

1.4576 

1.4443 

1.0091 

0.5885 

0.7 

0.6407 

1.3700 

1.2760 

1.0736 

0.7026 

o 

00 

0.7517 

1.2791 

1.1123 

1.1499 

0.8244 

0.9 

0.8711 

1.1873 

0.9584 

1.2388 

0.9552 

1.0 

1.0000 

1.0966 

0.8175 

1.3413 

1.0966 
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Table  21. 

Converged 

solution. 

B = tt/2,  MQA. 

t 

y 

z 

m 

u 

w 

0.0 

0.0000 

6.1911 

9.0703 

0.6825 

0.0000 

0.1 

0.0685 

6.1155 

8.8501 

0.6910 

0.1691 

0.2 

0.1387 

5.8977 

8.2309 

0.7165 

0.3424 

0.3 

0.2124 

5.5620 

7.3207 

0.7597 

0.5241 

0.4 

0.2913 

5.1422 

6.2573 

0.8217 

0.7188 

0.5 

0.377! 

4.6739 

5.1695 

0.9041 

0.9313 

0.6 

0.4729 

4.1888 

4.1521 

1.0088 

1.1668 

0.7 

0.5800 

3.7118 

3.2603 

1.1384 

1.4312 

0.8 

0.7015 

3.2600 

2.5149 

1.2962 

1.7310 

0.9 

0.8403 

2.8436 

1.9134 

1.4861 

2.0735 

1.0 

1.0000 

2.4674 

1.4406 

1.7126 

2.4674 
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Table  22.  Optimal  mass  distribution  m(t)/m(0)  for  several 
values  of  the  frequency  parameter  B . 


t 

II 

O 

3 = tt/6 

3 = 7t/5 

3 = 7t/4 

3 = 7t/3 

3 = u/2 

0.0 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

0.1 

1.0000 

0.9972 

0.9960 

0.9938 

0.9891 

0.9757 

0.2 

1.0000 

0.9891 

0.9843 

0.9757 

0.9573 

0.9074 

0.3 

1.0000 

0.9757 

0.9652 

0.9464 

0.9074 

0.8071 

0.4 

1.0000 

0.9573 

0.9394 

0.9074 

0.8431 

0.6898 

0.5 

1.0000 

0.9344 

0.9074 

0.8603 

0.7691 

0.5699 

0.6 

1.0000 

0.9074 

0.8703 

0.8071 

0.6898 

0.4577 

0.7 

1.0000 

0.8768 

0.8290 

0.7496 

0.6094 

0.3594 

0.8 

1.0000 

0.8431 

0.7845 

0.6898 

0.5312 

0.2772 

0.9 

1.0000 

0.8071 

0.7378 

0.6294 

0.4577 

0.2109 

1.0 

1.0000 

0.7691 

0.6898 

0.5699 

0.3904 

0.1588 
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Table  23.  Comparison  of  the  optimal  beam  (subscript  o) 
with  the  constant  section  beam  (subscript c)  . 


6 

I 

o 

I 

c 

I /I 
o c 

0 

0.00000 

0.00000 

1.0000 

0.0000 

7r/6 

0,30014 

0.30229 

0.9929 

0.0071 

n/5 

0.44954 

0.45650 

0.9848 

0.0152 

tt/4 

0.75458 

0.78539 

0.9608 

0.0392 

7r/3 

1.56091 

1.81379 

0.8606 

0.1394 

■it/2 

5.29597 

00 

0.0000 

1.0000 

i 
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Table  24.  Comparison  of  the  optimal  beam  (subscript  o) 
with  the  constant  section  beam  (subscript  c) . 


lo/Ic  ^^c’^o^/^c 


o 

c 

o'  c 

c o'  c 

0.0 

0.0000 

0.0000 

1.0000 

0.0000 

0.1 

0.0100 

0.0100 

1.0000 

0.0000 

0.2 

0.0405 

0.0405 

0.9999 

0.0001 

0.3 

0.0927 

0.0928 

0.9992 

0.0008 

0.4 

0.1687 

0.1691 

0.9976 

0.0024 

0.5 

0.2715 

0.2731 

0.9941 

0.0059 

0.6 

0.4053 

0.4104 

0.9874 

0.0126 

0.7 

0.5754 

0.5896 

0.9760 

0.0240 

00 

• 

o 

0.7887 

0.8237 

0,9575 

0.0425 

0.9 

1.0537 

1.1341 

0.9291 

0.0709 

1.0 

1.3810 

1.5574 

0.8868 

0.1132 

1.1 

1.7839 

2.1612 

0.8254 

0.1746 

1.2 

2.2784 

3.0865 

0.7382 

0.2618 

1.3 

2.8845 

4.6827 

0.6160 

0.3840 

1.4 

3.6263 

8.1170 

0.4468 

0.5532 

1.5 

4.5338 

21.1521 

0.2143 

0.7857 

Tr/2 

5.2959 

OO 

0.0000 

1.0000 
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10 . Conclusions 

In  this  report,  we  have  investigated  the  axial  vibra- 
tion of  a cantilever  beam  both  analytically  and  numerically. 

We  have  determined  the  mass  distribution  that  minimizes  the 
total  mass  for  a given  value  of  the  fundamental  frequency 
using  both  the  sequential  ordinary  gradient-restoration  al- 
gorithm (SOGRA)  and  the  modif ied-quasilinearization  algoritlim 
(MQA) . Concerning  the  minimum  value  of  the  mass,  SOGRA  leads 
to  a solution  precise  to  at  least  4 significant  digits  and 
MQA  leads  to  a solution  precise  to  at  least  6 significant 
digits. 

Comparison  of  the  optimal  beam  (a  variable-section  beam) 
with  a reference  beam  (a  constant-section  beam)  shows  that  the 
weight  reduction  depends  strongly  on  the  frequency  parameter  8. 
This  weight  reduction  is  negligible  for  6^0,  is  0.6%  for 
8=0.5,  is  11.3%  for  8=1,  is  55.3%  for  8=1.4,  and  approaches 
100%  for  8 -»■  Tr/2. 
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11 . Appendix;  Analytical  Solutions 

In  this  appendix,  we  summarize  the  analytical  solutions 
available  in  two  limiting  cases:  the  constant-section  beam 
and  the  optimal  beam. 

Constant-Section  Beam.  For  this  structure,  the  func- 
tions y(t),  z(t),  m(t),  u(t),  w(t)  are  given  by 


y = 

sin ( 3t) /sin3  , 

(91-1) 

z = 

3^cos  ( Bt) /cos3  , 

(91-2) 

m = 

3tan3  , 

(91-3) 

u = 

3cos ( 3t) /sin3  , 

(91-4) 

w = 

-B^sin ( 3t) /sin3  . 

(91-5) 

The  associated  mass  integral  is  given  by 

I = BtanB  . (92) 

Optimal  Beam.  For  this  structure,  the  functions  y(t), 
z(t),  m{t),  u(t),  w(t)  are  given  by  (Ref.  8) 

y = sinh(8t)/sinh(3)  , (93-1) 

z = e^cosh(3)/cosh(Bt)  , (93-2) 

m = 3sinh ( 3) cosh ( 3) /cosh^ ( 3t)  , 


(93-3) 
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u = Bcosh (8t) /sinh ( 3)  , (93-4) 

w = 3^sinh(3t)/sinh(3)  . (93-5) 

The  associated  mass  integral  is  given  by 

I = sinh^(3)  . (94) 

Limiting  Case.  For  3=0,  Eqs.  (91)  and  (93)  reduce  to 

y=t,  z = 0,  m=0,  u=l,  w=0,  (95) 

and  the  mass  integrals  (92)  and  (94)  yield 

1 = 0.  (96) 

Therefore,  for  3=0,  the  constant-section  beam  and  the 
optimal  beam  are  identical . 
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